Isostatic phase transition and instability in stiff granular materials. 
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Structural rigidity concepts are used to understand the ori- 
gin of instabilities in granular aggregates. It is first demon- 
strated that the contact network of a noncohesive granular 
aggregate becomes exactly isostatic when I = ke/ fi » 1, 
where k is stiffness, e is the typical interparticle gap and Jl is 
the typical stress induced by loads. Thus random packings of 
stiff particles are typically isostatic. Furthermore isostaticity 
is responsible for the anomalously large susceptibility to per- 
turbation observed in granular aggregates. The load-stress 
response function of granular piles is critical (power-law dis- 
tributed) in the isostatic limit, which means that slight over- 
loads will produce internal rearrangements. 



Photoelastic visualization experiments show 
clearly defined stress-concentration paths in non-cohesive 
granular materials under applied load. These often suffer 
sudden rearrangement on a global scale when the load 
conditions are slightly changed, evidencing a degree of 
susceptibility to perturbation not usually present in elas- 
tic materials. It is rather possible that this intrinsic in- 
stability be responsible for much of the interesting phe- 
nomenology of granular materials [3,4]. Recently a num- 
ber of phenomenological models [2,5-0 have been put 
forward, which succeed to reproduce several aspects of 
stress propagation in granular systems, and the issue of 
instability has been addressed by noting that the load- 
stress response function may take negative values H . It 
is the purpose of this letter to show that structural rigid- 
ity concepts help us understanding the origin of insta- 
bility in granular materials, linking it to the topological 
properties of the system's contact network. 

Structural rigidity || studies the conditions that a net- 
work of points connected by rotatable bars (representing 
central forces) has to fulfill in order to sustain applied 
loads. A network with too few bars is flexible, while if 
it has the minimum number required to be rigid it is 
isostatic. Networks with bars in excess of minimal rigid- 
ity are overconstrained, and are in general self-stressed. 
Concepts from structural rigidity were first introduced 
in the study of granular media by Guyon et al |To| , who 
stressed that granular systems are not entirely equivalent 
to linear elastic networks since in the former only com- 
pressive interparticle forces are possible. We next show 
that this constraint has far-reaching consequences for the 
static behavior of stiff granular aggregates. 

Consider a d-dimensional frictionless granular pile in 
equilibrium under the action of external forces Fi (gravi- 



tational, etc) on its particles. Imagine building an equiv- 
alent linear-elastic central-force network (the contact net- 
work), in which two sites are connected by a bond if and 
only if there is a nonzero compression force between the 
two corresponding particles. Because of linearity, stresses 
fij on the bonds of this equivalent system can be decom- 
posed as fij = f'j + f\° ad where /£• are self-stresses, 
and fl° ad are load-dependent stresses. These last are 
linear in the applied load and are do not change if all 
stiffnesses are rescaled. Self-stresses in turn do not de- 
pend on the applied load, but are linear combinations of 
terms of the form kijCij where kij are the stiffnesses of 
the bonds, and their length-mismatches. In a granu- 
lar pile, length-mismatches are due to interparticle gaps, 
and therefore will depend on the distribution of radii 
and on the characteristics of the packing. Furthermore, 
self-stresses can only arise within overconstrained sub- 
graphs P,|lO||, i.e. those with more contacts than strictly 
necessary to be rigid. It is easy to see that a bounded 
overconstrained subgraph with nonzero self-stresses must 
have at least one negative (traction) self-stress. It suf- 
fices to consider a joint belonging to the envelope of the 
overconstrained cluster: since bonds can only reach it 
from one side of the frontier, stresses of both signs are 
necessary in order for the joint to be equilibrated. Now 
rescale all stiffnesses according to k — > Afc. In doing 
so, self-stresses are rescaled by A, but load-dependent 
stresses remain constant. Thus if self-stresses were non- 
zero, in the limit A — > 00 at least one bond of the network 
would have negative total stress, which is a contradiction. 
Therefore stiff granular piles must must either: a) have 
zero length-mismatches, or b) have no overconstrained 
graphs at all. Condition a) cannot be satisfied if the par- 
ticles have random polydispersity, no matter how small, 
or if the packing is disordered. Therefore there can be no 
overconstrained subgraphs in polydisperse or disordered 
packings in the large-stiffness limit. In other words: the 
contact network of a granular packing becomes isostatic 
when the stiffness is so large that the typical self-stress, 
which is of order ke, would be much larger than the typ- 
ical load-induced stress Jl JH|/. 

The isostaticity condition above is perhaps simpler to 
understand when cast in the following terms: granular 
packings will only fail to be isostatic if the applied com- 
pressive forces are strong enough to close interparticle 
gaps establishing redundant contacts. 

Thus, real disordered packings will typically be iso- 
static (interparticle gaps are large) unless its particles 
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are strongly deformed by the load. This explains why 
the average coordination number of random sphere pack- 
ings is usually close to 6 (see and references therein) . 
Isostaticity was also reported in numerical simulations of 
rigid disks JFJ . 

Finally we note that an adimensional "isostatic- 
ity parameter" can be defined as J = ke/fi,, and 
that I >> 1 corresponds to the isostatic limit. 
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FIG. 1. Appropriately choosing among these three iso- 
static configurations for each site, only compressive stresses 
are produced on a triangular packing. First S is chosen with 
probability 1/2. If S is not chosen then either R or L are, 
depending on the sign of the horizontal force acting on the 
site. 

We now discuss the consequences of isostaticity for the 
static behavior of a pile. It is possible to obtain use- 
ful insight from recent studies of the related problem 
of central- force rigidity percolation [jig . Rigid back- 
bones are found to be composed of large overconstrained 
clusters, isostatically connected to each other by critical 
bonds (also called red bonds). Cutting one critical bond 
is enough to produce the collapse of the entire system, 
because each of these is by definition essential for rigid- 
ity. In percolation backbones though, the number of such 
critical bonds is not extensive, but scales at p c as L x ' v 
where v is the correlation- length exponent |r|. Thus, if 
we perturb (cut or stretch) a randomly chosen bond in a 
percolation backbone, most of the times the effect will be 
only be local since no critical bond will be hit. The new 
element in stiff granular contact-networks is the fact that 
all contacts are isostatic, or critic, i.e. there is extended 
isostaticity. Thus we may expect stiff granular systems 
to have a large susceptibility to perturbation since cut- 
ting (stretching) a bond will often produce a large part 
of the system to collapse (move). 

Let us now quantify these ideas. We perturb the sys- 
tem by introducing an infinitesimal change SI in the 
length of a randomly chosen bond, and record the induced 
displacement Sxi suffered by particle centers in equilib- 
rium. The system's susceptibility to perturbation is then 
defined as D = J2iLi \Sxi/Sl\ 2 . We propose to measure 
D(O v ) as a function of the density O v of overconstraints 
(excess contacts) randomly located on the network. Iso- 
static piles have O v = 0. 

A simple one-dimensional model for the propagation of 
perturbations can be analytically solved (Tt]] for arbitrary 
values of the density O v of overconstraints. For any non- 
zero Vl D as defined above takes a finite value, but 
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diverges as O v for O v — * 0. Therefore there is a phase 
transition at the isostatic point O v = 0. 



We now analyze a two dimensional system numeri- 
cally. In the spirit of previously studied models [H , we 
consider a triangular packing of H layers height, with 
one of its principal axis parallel to gravity, and made of 
disks with small random polydispersity SR and weight 
W. Since SR is small, disk centers are approximately lo- 
cated on the sites of a regular triangular lattice. If the 
stiffness k is large enough (/ = kSR/WH >> 1), the 
contact network will be isostatic. We enforce isostatic- 
ity in our model by letting each site be supported from 
below by only two out of its three neighbors. This gives 
three possible local configurations which are depicted in 
Fig.[|. By appropriately choosing among these |17|] , ran- 
dom isostatic networks with only compressive stresses are 
generated. In order to study the effect of a finite density 
O v of overconstraints (which would appear if the stiffness 
is lowered), we furthermore let all three bonds be present 
with probability O v at each site. 

After building a disordered network in this way, a ran- 
domly chosen bond in the lowest layer is stretched, and 
the induced displacement field is measured. After aver- 



aging over disorder, the stress distribution |17| is found 
to decay exponentially for large stresses, in accordance 
with previous work [p|Jl5|] . 

The results for the susceptibility D y (H,O v ) = 

Y^Li <%» 2 are shown in Fig. ||a. Here Syi is the vertical 
displacement of site i due to a unitary bond-stretching, 
as measured on packings of N = H x H particles. For 
O v > 0, D y goes to a finite limit for large sizes H, but 
diverges with system size if O v — 0. Measurements on 
isostatic packings of up to H = 2000 layers jl7| show that 
log 13 y oc H, i.e. the divergence of D y is exponential with 
size when O v = 0. Thus there is a surprising phase tran- 
sition at O v = 0, where anomalously large susceptibility 
sets in. 

In order to understand how displacements propagate 
upwards, we measure the probability distribution Ph(Sy) 
to have a vertical displacement Sy, h layers above the 
perturbation. Numerical results for isostatic systems 
with H = 2000 are shown in Fig. |^b. For large h, 
Ph(Sy) decays as a power-law with an /i-dependent cut- 
off: P h {5y) ~ h-P\Sy\- , Sy < S M (h). As seen in Fig. |b, 
Suih) grows exponentially with height h. This produces 
the observed exponential divergence of D y . Similar mea- 
surements were done on systems with a finite density of 
overconstraints O v , in which case the distribution of dis- 
placements presents a height- independent bound (lTj . 

The puzzling appearance of exponentially large dis- 
placements on isostatic piles can be explained as due to 
the existence of "lever configurations" or "pantographs" , 
which amplify displacements. Fig. || shows an exam- 
ple of a pantograph with amplification factor 2. Given 
that this and similar mechanisms appear with a finite 
density per layer, it is clear that the second moment of 
Ph(Sy) will grow exponentially with system height. Fur- 
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thermore, it is easy to understand why this amplification 
effect only exists in the isostatic limit: Pantographs as 
the one in Fig. || are no longer effective if blocked by over- 
constraints, for example if an additional bond is added 
between site A and the site below it. In this case, a 
stretching of bond B would induce stresses in the whole 
pantograph, and only a small displacement of site A. 
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tion H of the stretched bond. The network's total energy 
can be written as E = J2? =1 W lVt + 1/2 £ fe k b (l b - l° b ) 2 
where the first term is the potential energy and the sec- 
ond one is a sum over all bonds and accounts for the elas- 
tic energy. l b are bond lengths in equilibrium and I® their 
repose lengths. Upon infinitesimally stretching bond 
equilibrium requires that J2i ^i§i^ + J2 ov ^ovQov ~ 
lo V )*§ff = 0, where the second sum goes over bonds ov 
that belong to the same overconstrained graph as b' does. 
This is so since bonds not overconstrained with respect 
to b' do not change their lengths as a result of stretching 
b'. Since stress fb on bond b is /& = fc&(Z° — l b ) this may 
be rewritten as 



V f ^l = ^Tw ^- 



(i) 



20 40 60 80 100 



If b' does not belong to an overconstrained graph, 
the left hand sum only contains bond b' itself, there- 
fore fb — showing that, in the isostatic case 
(no overconstrained graphs at all), the induced dis- 
placement Sy^ = ^ is equal to the response func- 
tion of stress with respect to an overload on site i. 
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FIG. 2. a)Total susceptibility D y versus system height 
H in layers, as numerically measured on two-dimensional tri- 
angular packings. The fraction of overconstraints (fraction 
of sites supported by three lower neighbors) is: 0.00 (open 
circles), 0.01 (squares), 0.02 (diamonds), 0.05 (triangles) and 
1.00 (full circles), b) The probability Ph(5v) to have an in- 
duced vertical displacement 5y, H layers above the pertur- 
bation, as obtained in a numerically exact fashion for the 
isostatic (O v = 0) triangular piles described in the text. Re- 
sults are shown for H = 200, 400, 600, . . . , 2000. Only positive 
values of 8y have been plotted here. 

In order to formalize the relationship between these find- 
ings and the observed unstable behavior |l],|j of granular 
materials, we now demonstrate the equivalence between 
induced displacements and the load-stress response func- 




FIG. 3. The observed exponential growth of induced dis- 
placements is due to the existence of random "pantographs" 
as the one shown in this figure. Upon stretching bond B by 
an amount 5, site A moves vertically by an amount 28. Con- 
versely a unitary weight at A produces a stress of magnitude 
2 on bond B. This is a consequence of the general equivalence 
between induced displacements and the load-stress response 
function. 

Taking averages with respect to disorder we obtain < 
fb >= J2i Wi < Sy^ > and since average stresses grow 
as H we must have < Sy^ >~ if -1 . We thus see that 
there must be delicate cancellations in P(5y), since its 
second moment diverges as exp{7J} while its first mo- 
ment goes to zero with H. This shows that dy (and 
therefore the response function) takes exponentially large 
values of both signs, (P(Sy) is approximately symmetric). 
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Thus a positive overload at site i would often produce a 
(very large) negative stress on bond b, implying the need 
for rearrangement since negative stresses are not allowed. 

The existence of negative values for the response func- 
tion was first discussed in relation with instability, in 
the context of a phenomenological vectorial model for 
stress propagation f|. The results of the present work 
demonstrate that the response function takes exponen- 
tially large negative values, and the system is unstable, 
because of the isostatic character of the contact network. 

To conclude, we have shown that granular packings are 
exactly isostatic when / = fee/ fi, is much larger than one, 
which holds for typical disordered packings, and also for 
stiff enough regular packings with random polydispersity. 

For isostatic packings, the distribution of displace- 
ments induced by a perturbation is power-law with an 
exponentially large cutoff. A susceptibility to pertur- 
bation can be defined, which diverges upon increasing 
I. Thus, an isostatic phase transition takes place in the 
limit of large /. 

Induced displacements were furthermore shown to be 
equivalent to the load-stress response function of the per- 
turbed bond. Our results for induced displacements thus 
mean that response functions take exponentially large 
values, as well positive as negative, in the isostatic limit. 
This explains why stiff granular piles are unstable. Any 
non-zero density of overconstraints destroys criticality 
and therefore instabilities will not be present when the 
isostaticity parameter / is small. I can be reduced by 
reducing the stiffness, the interparticle gaps, or by in- 
creasing the load. 
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